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ABSTRACT 

We present a gravitational lcnsing analysis of the galaxy cluster Abell 1689, incorporating measure- 
ments of the weak shear, flexion, and strong lensing induced in background galaxies. This is the first 
time that a shapclet technique has been used to reconstruct the distribution of mass in this cluster, 
and the first time that a flexion signal has been measured using cluster members as lenses. From weak 
shear measurements alone, we generate a non-parametric mass reconstruction, which shows significant 
substructure corresponding to groups of galaxies within the cluster. Additionally, our galaxy-galaxy 
flexion signal demonstrates that the cluster galaxies can be well-fit by a singular isothermal sphere 
model with a characteristic velocity dispersion of a v = 295 ± 40 km/ s. We identify a major, distinct 
dark matter clump, offset by 40ft.~ 1 kpc from the central cluster members, which was not apparent from 
shear measurements alone. This secondary clump is present in a parametric mass reconstruction using 
flexion data alone, and its existence is suggested in a non-parametric reconstruction of the cluster using 
a combination of strong and weak lensing. As found in previous studies, the mass profile obtained by 
combining weak and strong lensing data shows a much steeper profile than that obtained from only weak 
lensing data. 

Subject headings: cosmology: observations - galaxies:clusters:general - galaxies:photometry - 
gravitational lensing 

1. INTRODUCTION 

Gravitational lensing studies provide a powerful tool for mapping out the surface mass density of galaxies and clusters 
of galaxies (e.g. Bradac et al. 2005ab, 2006; Natarajan & Springel 2004; Diego et al. 2005a; Cacciato et al. 2006; 
Abdelsalam et al., 1998). In particular, cluster lensing studies provide constraints on both structure formation models 
and the mean mass density of the universe. Additionally, accurate knowledge of the mass distribution in clusters enables 
better determination of the relationship between dark matter, gas and galaxies within the cluster, and allows us to probe 
how well light traces mass. For example, there has been tremendous excitement recently about maps from gravitational 
lensing of the distribution of baryonic mass and dark matter in large-scale structure (Massey et al. 2007a) and particularly 
the cluster 1E0657-56 (the "Bullet Cluster"; Bradac et al. 2006), in which there is a marked offset between the peaks of 
baryonic and non-baryonic mass. 

Traditional approaches to weak gravitational lensing focus on estimating the mass of the lens by measuring the induced 
cllipticity (shear) in images of background galaxies (see, e.g. Kaiser, Squires & Broadhurst, 1995 [hereafter KSB], for 
a review of the canonical approach inverting cllipticitics; Bartclmann & Schneider, 2001, provide an excellent review of 
the subject). Underlying these approaches is the assumption that any variations in the lensing field over the scale of the 
galaxy image may be neglected. Flexion has recently been introduced as a means to describe higher order lensing effects 
which arise as a result of small-scale variations in the lensing field (Goldberg and Bacon, 2005 [hereafter GB05]; Bacon, 
Goldberg, Rowe and Taylor, 2006 [hereafter BGRT], Goldberg & Leonard, 2007 [hereafter GL07]). The flexion technique 
enables us to probe structure on smaller scales than does a shear analysis, and is particularly useful in galaxy- galaxy 
lensing studies. 

In this paper we present several mass reconstructions of Abell 1689, a rich cluster of galaxies at a redshift of ~ 0.18, 
using images taken by the HST Advanced Camera for Surveys (ACS). This cluster has been well studied in the context 
of both weak and strong lensing, and various parametric models for the mass distribution have been tested (see, for 
example, Bardeau et al., 2005; Halkola, Seitz & Pannella, 2006; Saha, Read & Williams, 2006; Zekser et al., 2006; 
Umctsu, Broadhurst, Takado & Kong, 2005; Broadhurst et al., 2005a, Diego et al., 2005b, Sharon et al, 2005; King et 
al., 2002; King, Clowe & Schneider, 2002; Dahlc, Kaiser, Squires & Broadhurst, 2001; Dye et al., 2001; Taylor & Dye, 
1998; Tyson & Fischer 1995). Recently, Halkola, Seitz & Pannella (2007) have investigated the size of individual galaxy 
haloes within the cluster using strong lensing, which is typically difficult to do in clusters using traditional weak lensing 
techniques. 

Likewise, a number of researchers, including the original observing team (Broadhurst et al. 2005b; see also Diego et 
al. 2005b) have produced parametric and non-parametric mass estimates using strong lensing from the same ACS images 
in our sample. Most recently, Limousin et al. (2006 [hereafter L06]) have used these images, as well as images from the 
CFH12k camera, to model the cluster using both strong and weak lensing measurements, and find good agreement in the 
mass estimates derived from these two sets of measurements. 

We present here a non-parametric mass reconstruction of this cluster from shear measurements, as well as a galaxy- 
galaxy lensing study from measurements of flexion. It is typically very difficult to extract a weak shear galaxy-galaxy 
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signal in clusters of galaxies, as the smooth component of the cluster field dominates. However, flexion is better suited 
to detect perturbations to this lensing field on smaller scales, and thus is particularly adept at picking out lensing signals 
from individual galaxies in the cluster. 

We show that our flexion data can be used to construct a parametric mass reconstruction that successfully identifies 
substructure within the cluster. Additionally, we combine our non-parametric shear mass reconstruction with a non- 
parametric strong lensing reconstruction. The combined convergence map shows significant substructure, and suggests 
the presence of the secondary dark matter clump described by L06. Additionally, this convergence map shows a steeper 
radial profile in the center of the cluster than the one generated by shear measurements alone. 

This paper is structured as follows: we begin in § [2] with a brief review of flexion formalism and a discussion of our 
measurement techniques, including an overview of the formalism underlying the shapelets and HOLICs methods. We 
describe our data and processing pipeline in detail in § [3] and present the results of our study in § [4] 

2. MEASURING SHEAR AND FLEXION 
2.1. Flexion Formalism 

It is helpful to begin by specifying our notation. Consider an image (at z = oo) which, in the absence of lensing would 
be observed at angular position (3. This may be related to the observed coordinate, 6 via the 2nd order transformation: 

A — AijOj + -DijkOjOk, (1) 
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These operators give rise to an asymmetrical distortion in the image, namely a skewness and a bending or arciness, as 
well as a shift in the centroid of the light distribution. BGRT define two "flexions" as follows: 



T = 0*7 = (0i 7i + 0272) + «(0i72 - 027i) (5) 

and 

Q = 07 = (0i 7i - 0272) + «(0i72 + 027i), (6) 

where d — d\ + i&2 is the complex gradient operator defined in BGRT. First flexion, J 7 , has an m — 1 rotational 
symmetry (i.e. it is a vector) and gives rise to a skewness in the light distribution of a galaxy image, as well as a shift 
in its centroid. It is also simply related to the convergence, k, by T — On. Second flexion, Q, has an m — 3 rotational 
symmetry, and gives rise to bending in the lensed image. 

2.2. Measurement Techniques 

Two distinct methods have been proposed for measuring flexion in real images. The method described in GB05 and 
BGRT involves the decomposition of galaxy images into shapelets (Rcfrcgicr 2003), followed by an "active" perturbation 
of coefficients (adopting the terminology of the Shear TEsting Programme; Heymans et al., 2006; Massey et al., 2006). 
A "passive" method, involving high-order moments of a galaxy's shape, referred to as Higher Order Lensing Image's 
Characteristics (HOLICs), was first related to flexion by Okura et al. (2007). These techniques have been refined, 
compared, and discussed in detail in GL07. We also note that others (Irwin & Shmakova 2005; 2006) have proposed a 
statistic similar to flexion, using a x 2 minimization technique to fit a polynomial radial profile with first- and second-order 
perturbations. 



3 



2.2.1. Shapelets 

Several researchers (Refregier, 2003; Bernstein & Jarvis, 2002) have proposed techniques for decomposing images into 
"shapelets", simple basis functions composed of the two-dimensional, Gaussian- weighted Hermite polynomials B nm (6). 
Any isolated image f(0) can be expressed as a sum 

f(0) = ^2f nm B nm (G) , (7) 

mn 

where the weights f nm are known as "shapelet coefficients". As shown in Refregier (2003), Refregier & Bacon (2003), 
Massey & Refregier (2005) and Masscy ct al. (2007b), shapelets are an especially useful basis for our purposes because 
image transformation operations induced by gravitational lcnsing typically produce very compact transfers of power 
between shapelet coefficients. 

In the shapelets framework, the lensing operators can be expressed elegantly in terms of the quantum mechanical raising 
and lowering operators, a and a\ which are combinations of the 0i and di operators. 
We can express the lens equation as: 

f={l + Kk + lj S j + li>j S ( g ) )f, (8) 

where /' represents the unlensed image, / refers to the lensed image, and the various operators K, Sj and Sh' are 
expressed in terms of a and (the reader is referred to GB05 for the details) . 

An important feature of the operators in Equation 8 is that the K and Sj operators transfer power in shapelet space 

— * (2) 

such that |An| + \ Am\ = 2, whereas the second-order operators $y yield |An| + |Am| = 1 or 3. This provides a 
straightforward mechanism for extracting the first and second order signals. 

As with a shear analysis, our flexion analysis makes the assumption that any intrinsic flexion signal will be randomly 
oriented, and thus will average to zero. This means that we expect that on average the (n + m =) odd shapelet modes 
will have a zero signal. 

In practice, the shapelet series will need to be truncated, in which case there are two parameters that must be specified 
for the decomposition. These parameters are n max , the maximum order of the shapelet series, and (3, the characteristic 
scale of the shapelets: 

ZW0)«exp(~^_if). (9) 
As discussed in Refregier (2003), the optimal choices for f3 and n max are: 

— V ^ @max @min 
9 max -. 

n m ax ^ p 1, (10) 

where 9 max and 9 m i n are the maximum and minimum scales on which one expects to resolve structure in the image, 
respectively. For the data used here, we find 9 min — 1.0 pixels and 8 ma x — 1.5-\/a 2 + b 2 give good shapelet reconstructions 
of our galaxy images, where a and b are the semi-major and minor axes of the object as measured during source extraction, 
respectively. 

Images also need to be corrected for PSF and pixellisation effects, which both tend to dilute the lensing signal. One 
of the advantages of the shapelets technique is that an explicit PSF deconvolution, and an integration within pixels, can 
be carried out prior to any parameter estimation. Refregier & Bacon (2003) give an explicit form for the deconvolution, 
and Massey & Refregier (2005) derive the form of the integrals. We perform these tasks using the publicly available IDL 
software from the shapelets web pag^j] 

2.2.2. HOLICs 

We have discussed above how shapelets may be used to measure both the shear and the flexion of a lensed image. 
Historically, however, shear has primarily been measured using combinations of weighted moments of the image. Likewise, 
the HOLICs technique proposed by Okura et al. (2007), and subsequently refined in GL07, provides a straightforward 
and fast way of estimating flexion using weighted moments of the light distribution of galaxies. Okura et al. define the 
following complex terms: 

^ _ (Qui + Q122) + i(Qii2 + Q222) ^ 

and 

, (Qui - Q122) + «(3Qii2 - Q222) n n , 

6 = ^ , ( 12 J 

1 http: / / www. astro, caltech .edu / ~ rj m/ shapelets / 
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where 

£ = Qllll + 2Qii22 + Q2222- (13) 

As pointed out in GL07, the relationship between the above quantities and flexion estimators is best expressed as a matrix 
equation: 



M 
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(14) 



where M. is a 4 x 4 matrix whose elements are proportional to linear combinations of Qiju and QijQu- The reader is 
referred to GL07 for the explicit forms of these elements. 

Equation [14] describes how flexion could be estimated from unweighted moments in the absence of any measurement 
noise. However, when dealing with real images, and particularly those dominated by sky noise, measurement of these 
moments is inherently noisy. It is thus necessary to include a weighting function in our measurements of the moments, 
which gives higher weighting to the central region of the postage stamp of a galaxy (where the actual image lies) than to 
the extreme regions (which we expect to be dominated by noise). 

We use a Gaussian filter: 
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and define the weighted moments as: 

Qij = 4 f(0i - 0l)(9j - 9~)f(6)W(9)d 2 e. (16) 

Using the weighted moments necessitates two corrections to Equation [14] The first of these corrections results from 
the fact that there is a difference in the centroid shift between the weighted and unweighted calculations. The second 
correction has to do with the fact that the total flux is not conserved by lensing. Ordinarily, this is related by the Jacobian 
of the coordinate transformation, however when using a window function, the transformation must be considered explicitly. 

Incorporating these corrections, we can write 

M = M(Q ij ,...) + AM, (17) 

where AM is a 4 x 4 matrix representing the correction terms. These terms are proportional to sums of Qijkimn and 
QijQkimn- The reader is referred to GL07 for details regarding the computation of these terms, which are expressed in 
full in the Appendix to that paper, and available as an IDL function from the flexion websit^] 

3. DATA AND PROCESSING PIPELINE 
3.1. Data 

Our data consists of 20 HST ACS images of Abell 1689, taken using the Wide Field Camera (WFC) during HST cycle 
11 by H. Ford. These images are all 2300-2400 second exposures covering a square field of view of 3.4' on each side 
(corresponding to ~ 460/i _1 kpc at the distance of the cluster), and are described in detail in Broadhurst et al. (2005b). 
Of these images, 4 were taken using the F475W filter (G-band), 4 using the F625W filter (R-band), 5 using the F775W 
filter (I-band) and 7 using the F850LP filter (Z-band). 

Due to the sensitivity of the ACS camera, these exposures resolve objects down to a magnitude of ~ 27 in the ri- 
band, which makes them ideal for looking at very faint background sources. Additionally, this sensitivity gives rise 
to a very high background source count: our final catalog had ~ 200 sources/ arcmin 2 . However, in images of Abell 
1689, the central cluster galaxies are a dominant feature, with R-band magnitudes as low as 15-16. Figure [l] shows the 
magnitude distribution for sources in an R-band image, and two distinct peaks can be identified. The lower magnitude 
peak corresponds primarily to central cluster members and bright stars, while the background galaxies are seen in the 
larger, high-magnitude peak. 



2 http: //www. physics. drexel.edu/~goldbcrg/flexion/ 
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Fig. 1. — The figure shows the distribution in magnitudes of galaxies in an R-band image of Abell 1689. There are two 
clearly identifiable peaks in this distribution. While the majority of the objects in our sample are faint, with a magnitude 
peak around 26, there is also a low magnitude peak at around 19. This corresponds primarily to cluster members, which 
are concentrated towards the center of the images, and bright stars in the field 



The large spread in magnitudes seen in Figure [T] presents several difficulties in the source detection and extraction 
process. Typical source extraction software packages, such as SExtractor (Bertin & Arnouts, 1996), require various 
detection and object size thresholds to be set prior source extraction. In images where there is such a large range of 
brightnesses, setting lenient detection criteria will result in detection of fainter background sources, but will also lead to 
excessive blending of images in the neighbourhood of the brightest sources. On the other hand, a strict set of criteria 
will reduce the amount of blending seen, but will also favor detection of foreground objects over the background objects 



that we are interested in. In § 3.2.2 we discuss the implementation of a two-pass source detection strategy carried out 
using SExtractor. This technique is modelled on the "hot" and "cold" source detection strategy described by Rix et al. 
(2004), who first noted that a single set of SExtractor parameters is generally insufficient for the detection of all objects 
of interest in ACS images. 

3.2. Catalog Generation 
3.2.1. Co- addition of Images 

For the initial source detection, we co-added our images using the SWarp software packag^] This software produces 
background-subtracted, median-stacked images, which allows us to clean our images of spurious hotspots and bad pixels. 
Images are read in by the software, and the background is estimated and subtracted out. The images are then resamplcd 
and projected into subsections of the output frame using any of a number of astrometric projections defined in the WCS 
standard. Of the projections included in the software, we opted to use the distorted tangential (gnomonic) projection. 
This projection introduces very little distortion in images smaller than 10 degrees (see Goldberg & Gott, 2006), and is 
recommended in the SWarp software manual for small fields. The images are then co-added by taking the median value 
at any given location. 

In addition to cleaning our images of bad pixels, stacking the images increases the signal-to-noise within a given image, 
thus better enabling us to detect faint background sources. A single stacked image was generated in the G-, R- and 
I-bands. In the Z-band, we opted to create two stacked images. Our initial dataset consisted of 7 exposures in this band. 
For the purposes of source detection, we found that very little advantage was gained by stacking 7 images compared to 
stacking 3 or 4 frames. Thus, we created two Z-band images, comparisons between which provide an important test of 
the effectiveness of our source detection and extraction strategy. 

3.2.2. SExtractor Two-Pass strategy 

We use SExtractor (Bertin & Arnouts, 1996) to identify and extract sources in our images. This software allows the 
user to specify a number of input and output parameters, most notably a detection threshold and a minimum number 
of connected pixels required to be considered an object. Additionally, the software is able to produce background and 
variance maps for the input images, and to generate background-subtracted frames. 

As noted above, a single SExtractor run will generally result in an incomplete catalog of background sources. Thus, 
source extraction is carried out in two stages. The first stage is designed to detect only foreground objects by employing 
the cross-identification utility included in the software. This utility allows the user to supply a catalog of object positions, 
and to instruct the software either to include or exclude the sources in the catalog. We supplied the locations of known 
foreground objects compiled from the list of spectroscopically confirmed cluster members presented by Due et al. (2002), 



3 http: / /terapix. iap.fr/soft/swarp 
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prominent stars identified by eye, and objects identified as cluster members by the NASA/IPAC Extragalactic Database 
(NED). For this pass, we set the detection threshold at 4cr, and the minimum area of detected objects at 10 pixels. 

During this first run, background and variance maps were generated, and this information was used to mask out the 
foreground objects in order to simulate an emptier field. A friends-of-friends algorithm was used to identify connected 
pixels above our detection threshold at the locations of the foreground objects. These pixels were then masked out 
according to f pixel — /background + TZ-Vpixei where 1Z is a random number drawn from a standard normal distribution 
(jH = 0, a = 1), <7 p i xe i is the standard deviation of the background measured at the pixel location, and /background is 
the background level. Thus, the foreground objects are replaced with a simulated noisy background which has the same 
statistical properties as the background measured in that region of the image. 

A secondary SExtractor run was then performed on the masked image, using less stringent detection criteria (detection 
threshold: 2a, minimum area: 15 pixels), to pick out background galaxies in the field. This strategy alleviated any 
blending problems we might otherwise have had by removing foreground objects prior to lowering the detection threshold. 
Additionally, it removed from our catalog known foreground objects which, if included in the subsequent analysis, could 
dilute the measured shear and flexion signals. 

A final catalog of background objects was generated by comparing the image detections across different bands. An 
object was included in the catalog if it was detected in at least two different bands, one of which was either the R- oi- 
l-band. 

Whilst the stacked images allow us to better identify legitimate background sources, they are generally not suitable for 
carrying out measurements, particularly of shear. This is because the PSF varies from image to image; each exposure 
is offset slightly from the others, and this results in a very complicated PSF in the stacked image that is impossible to 
model simply. For this reason, the shapelet measurements were carried out on the unstacked images, allowing an explicit 
PSF deconvolution to be carried out during the shapelet analysis. 

It was thus necessary to run SExtractor on the individual frames, since the shapelet parameters used in the analysis 
depend on SExtractor measurements of galaxy shapes in a given image (9 max = 1.5y/a 2 + b 2 ). We used the two-pass 
strategy described above, using detection thresholds of la less than those used on the stacked images to reflect the lower 
signal-to-noise in the individual frames (the minimum source areas remained the same). Additionally, in the second 
pass, we required that SExtractor detect only those objects that were detected in the stacked frames and included in the 
master background object catalog. This requirement avoided any spurious detections resulting from bad pixels within the 
individual exposures. 

3.3. Shapelet Analysis and Correction of Image Distortions 

Source extraction was carried out on each individual frame, generating a catalog for that frame that contained infor- 
mation about the size, shape and magnitude of each detected object. The next task was to generate a postage stamp of 
each individual galaxy, and perform a shapelet decomposition in order to measure the shear and flexion. 

To generate a postage stamp, a circular region around each galaxy image was extracted. The size of this region was 
determined by the SExtractor estimate of the size of the galaxy. A segmentation map was then created by identifying all 
pixels belonging to the primary source, as well as any other sources found within the extracted region. At this stage, if the 
primary source was found to extend close to, or beyond, the boundary of the postage stamp, the size of the extracted region 
was increased appropriately. The segmentation map was then checked for blending of the primary source. Where there 
was no blending, any other sources in the region were masked out, using a similar procedure to that used in masking out 
the foreground objects. Postage stamps in which blending of the primary source was found were automatically rejected. 

The next step in the process is a shapelet decomposition and PSF deconvolution. As discussed in GL07, the computing 
time required for a shapelet decomposition scales roughly as 0^ aa . (not including a PSF deconvolution, which slows the 
process further). This means that for very large objects, the decomposition time becomes prohibitive. GL07 discuss a 
method for speeding up the decomposition in larger objects by re-gridding them into larger pixels. This method works 
adequately, however it does not take into consideration any PSF effects. We discuss the effect of the PSF on shear 
and flexion measurements below, and describe our strategy for minimizing the shapelet analysis time without significant 
reduction in the accuracy of our shape measurements. 

3.3.1. PSF Modelling Using TinyTim 

We modelled the PSF using the TinyTim software package (Krist, 1993). This software was designed to simulate the 
PSF of the WFC using a three-step process that allows it to account for variations in the PSF due to chip location and 
filter wavelength. We used the default WFC settings for the focus and PSF size in our TinyTim simulations. 

The ACS wide field camera has a well-studied geometric distortion in its images, which is a result of the off-axis location 
of the camera on the Hubble Space Telescope. This distortion has been modelled (Meurer et al., 2003), and is well fit 
by a 4 th order polynomial. TinyTim attempts to take this distortion into account by applying it to the PSF in the third 
stage of the PSF generation. 

This is necessary when correcting for the PSF at the point of parameter estimation. However, our shapelet technique 
involves an explicit PSF deconvolution before any shape measurements are made. Including this correction to the PSF 
prior to deconvolution does not amount to correcting the image for the geometric distortion. Thus we did not include this 
step in our PSF generation; rather, we used an undistorted PSF model in the deconvolution. At the point of parameter 
estimation, we applied a correction for the geometric distortion as described in § |3.4| 
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As the decomposition time in larger objects is of concern, one might ask what effect this model PSF has on shear and 



flexion measurements. The flexion measured in the PSF is small (typically of order 10 
a flexion in the image which scales roughly as: 



10 5 /pixel). This will induce 
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where apsF and a source refer to the semi-major axis of the PSF and the source, respectively. This scaling relation 
follows from equations 43-45 in GL07, and is discussed in more detail in Appendix [A] Clearly, this drops off very rapidly 
with increasing source size, thus the PSF will be a sub-dominant effect in measurements of flexion for images with a 
semi-major axis comparable to that of the PSF, which, for a typical R-band simulation, is found to be o>psf — 2.0 pixels. 
Thus, the flexion measuremen ts, in general, do not need to be corrected for PSF effects, particularly if a minimum size 
criterion is introduced (as in § 3.7). 



In order to assess how the PSF affects the measured ellipticities of galaxies, we simulated gaussian "galaxies" of various 
sizes and ellipticities, convolved them with a model R-band PSF in real space, and computed the change in ellipticity 
of the source. The ellipticity was computed using the unweighted moments of the light distribution, and we defined the 
change in ellipticity of the source as: 

A|e| S) /(4 0) -ei c) ) a + (4 0) -4 e) ) a (19) 
where the superscript (0) and (c) refer to the unconvolved and convolved images, respectively. 

The left panel of Figure jij shows the fractional change in ellipticity, , plotted as a function of n max (defined in 
2.2.1 ) for an elliptical source with | e |= 0.2. We also show in Figure [2] a comparison of the absolute change in ellipticity 



for a circular and an elliptical source. 



Fractional Change in Ellipticity as a Function of n_max 
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Fig. 2. — The left panel shows the fractional error in the ellipticity induced by the PSF plotted against n max = e ma * — 1 
for a series of test galaxies with an elliptical gaussian profile with | e |= 0.2. The right panel compares the absolute change 
in ellipticity for a circular source (dashed line) and an elliptical source with | e |= 0.2 (solid line). 

From Figure [2] we can see that the effect of the PSF on the ellipticity of the source drops to below 10% at n max = 20. 
Thus we defined a "large object cutoff" : any object with an n max > 20 was considered to be sufficiently large that 
PSF effects do not dominate, and an explicit deconvolution was not carried out. Additionally, as described in GL07, 
images with n max > 50 were re-gridded into larger pixels of size = n max /A0 in pixel units, in order to further reduce the 
decomposition time for very large objects. In order to exclude objects for which the shapelet fit was not successful, we 
rejected any sources for which the decomposition had a \ 2 per degree of freedom that was greater than 2.0. 

3.4. ACS Geometric Distortion Correction 

As described above, the ACS wide field camera induces a geometric distortion in images that is well fit by a 4 th order 
polynomial (Meurer et al., 2003): 

4 i 

i=0 j=0 
4 i 

i=0 j=0 
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where the primed coordinates refer to the undistorted frame, and the unprimed coordinates are the distorted pixel 
coordinates relative to the central pixel in each chip, aij and bij are the best fit coefficients, and are specific to each of 
the two chips on the ACS camera. The fit parameters can be found in the TinyTim software package. 

We can use these transformations to evaluate the A and D operators defined in Equations [2] and [3] This allows us to 
evaluate the shear and flexion induced by the geometric distortions as a function of location on the chip. These induced 
distortions can then be subtracted from the measured shear and flexion. 

We found that the mean magnitude of the induced ellipticity across the field of view was significant, with a mean of 
| 7 | = 0.0405, and a standard deviation of <7| 7 | = 0.0083. The induced flexion, however, was found to be negligible over 

the entire field of view, with \ Y\ = 3.05 x 10~ 4 , \ Q\ = 1.00 x 10~ 4 , a\ F \ = 4.61 x 10~ 5 and a\g\ = 2.34 x 10" 5 . 

In light of this, only the measured ellipticity was corrected for the geometric distortions. Additionally, since the PSF 
was found to have little effect on the flexion for images with a ~ 2.0, we opted to measure the flexion on the stacked 
images, for which the signal-to-noise is higher, and incorporated a lower limit on the size of images analyzed. We note 
that, as a result of the re-projection of the images during the stacking process, the geometric distortions are reduced in 
the stacked frames. 

3.5. Non-parametric Convergence Map Generation 

Each background galaxy had from one to twenty independent ellipticity estimates, depending upon how many of fields it 
was detected in. For those detected at least twice, we estimated the mean and measurement error of the galaxy ellipticities 
via the relation: 

_ {{ei-e)Y' 2 (m 

The mean and error were actually measured iteratively, and all measurements outside of 2.5cr of the mean were discarded 
as outliers. An ellipticity error was then assigned to each galaxy as a quadratic sum of the measurement uncertainty and 
0.3, the intrinsic standard deviation in ellipticities. 

We then binned the field into cells approximately 20 arcseconds on a side with an estimated ellipticity given by 
the weighted mean of the measurements in the individual frames, and with bin errors estimated using standard error 
propagation. Various rejection criteria in the pipeline (most notably blended sources and those with poor shapelet 
reconstructions or fewer than two detections) resulted in a reduced source count of 88/arcmin 2 , which is much higher 
than that typically found in weak lensing studies of this cluster (sec, for example, Broadhurst et al., 2005a). 

We then applied the finite inversion technique advocated by Seitz & Schneider (1995). They show that the density field 
may be solved iteratively (up to the mass sheet degeneracy) as: 

k{0) - k q = i J d 2 6[l - K{0')]n[V*{0 - e'){e){6% (22) 



where the convolution function T> can be written as: 



l-(l+|)exp(|) 

ie 2 y 



T>m=-- \n 'L , 2 — . ( 23 ) 



and our smoothing scale, 9 S , was set to half a binsize. 

The above method produces a reasonable mass reconstruction (see § |4| , however a better-constrained mass model can 
be derived from combining both weak and strong lensing methods (sec, for example, Kncib ct al., 2003; Smith et al, 2005; 
Bradac et al., 2005ab; 2006). In order to do this, we follow the method of Bradac et al. (2005a). Our multiple image 
catalog is taken from L06, and we assume that Di s /D s — 0.625, i.e. the mean redshift of the background sources is 0.9, 
as in L06. Note that applying the Bradac et al. technique with only weak lensing constraints produces a very similar 
reconstruction to that found using the approach of Seitz & Schneider (1995). 

3.6. Parametric Convergence Map Generation 

Flexion typically dominates the lensing signal on a much smaller scale than shear. Thus, the binning scale described 
above is not appropriate for flexion measurements as the signal would be entirely dominated by shot-noise. In this context, 
discrete estimates of the flexion of individual galaxies (parametric reconstructions) are in order. 

As the flexion was measured on the stacked images, each source had up to five flexion measurements. We included only 
those objects measured in at least 3 of the 5 stacked frames, and combined the flexion estimates using a similar iterative 
statistic to that described above for shear. Additionally, we only included objects that had a semi-major axis larger than 
0.12" (2.4 pixels), which were in the brightest 90% of the sample of background objects (which still yielded approximately 
75 galaxies/arcmin 2 ), and which had absolute values of flexion that were smaller than 4x the random scatter, to exclude 
extreme outliers. 

We found that the 2nd flexion produced significantly higher scatter than the 1st flexion, and thus, 1st flexion was the 
only one used in our analysis. 

To generate a convergence map, we modeled each cluster galaxy under the assumption that it is a singular isothermal 
sphere with the flexion profile given in BGRT: 
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F eW) = ■ ( 24 ) 

The flexion data was used to fit a% for each foreground galaxy for which at least three background sources were identified 
with an angular separation of < 10". Out of our initial sample of 62 foreground galaxies, only 36 met this criterion, and 
thus our final mass reconstruction was based on the fits obtained for these 36 galaxies. 

Since this is a noisy measure, some estimators predicted a negative value of <r^. These were retained in our map 
in order to keep the convergence estimate unbiased. The mean velocity disperson computed this way was found to be 
< a v >= 321fcm/s. After estimating the velocity dispersions, a density field was layed down: 

and then smoothed using a Gaussian kernel on a scale of 10". 

3.7. Galaxy- Galaxy Flexion Measurements 

We also used our flexion measurements in a galaxy-galaxy lensing study. Galaxy-galaxy measurements of shear are 
extremely difficult in cluster environments since a single source may be partially lensed by many different foreground 
objects. Additionally, the shear signal tends to be dominated by the smooth component of the cluster mass distribution, 
rather than the substructure. The flexion signal, however, drops off much more quickly with separation, and thus will 
almost always be produced by a single source. 

We carried out a pairwise comparison between each of the background images and the 62 cluster members in our 
foreground object catalog. For each, we considered the relative orientation angle and estimated "B" and "E" field flexion 
signals via the relationships: 

Te = T± cos(cf>) + T 2 sin(^) 

J~b = —T\ sin(0) + T 2 cos(4>), (26) 

where is the position angle of the background source with respect to the foreground galaxy. Lensing naturally gives rise 
to Te < and Tb — 0, thus the scatter in the measured B-modc signal gives an estimate of the noise in the measurements. 
For each, we estimated an uncertainty in the flexion via the relation: 



0.029\ 2 



j +°La. (27) 



where the former term on the right represents the scatter in intrinsic flexion (see GL07), and the latter represents the 
measurement error as determined from the scatter between frames. Weighting each data point by its signal-to-noise, we 
computed the average flexion signal as a function of lens-source separation. 

Finally, we performed least squares fitting to an assumed isothermal sphere profile as described in Equation [24] Other 
fits could be done, of course, but the isothermal model was both simple and seemed to fit the data well. 

Note that, in this case, we are considering the average signal over all foreground-background pairs, rather than fitting 
each individual foreground galaxy. 

4. RESULTS AND DISCUSSION 
4.1. Weak Shear Mass Reconstruction (Non-Parametric) 

Most mass reconstructions of Abell 1689 have either been parametric (e.g. L06; Halkola, Seitz & Panella, 2006; 
Broadhurst et al., 2005b; King, Clowe, & Schneider, 2002), have measured circularly averaged shears (e.g. L06; Bardeau 
et al., 2006; Broadhurst et al., 2005a; Umetsu et al., 2005), or have been non-parametric reconstructions using strong- 
lensing data (e.g. Diego et al., 2005b). Thus, it is interesting to compare weak non-parametric reconstructions to the 
above methods to verify that they produce consistent results. This would also lend fuel to the recent efforts to unify weak 
and strong lensing studies in order to further reduce uncertainty. 

In order to generate radially averaged shear profiles of the cluster, the authors mentioned above have made use of 
ground-based data, either exclusively or as a supplement to the ACS images of the cluster. Part of the reason for this is 
that the core of the cluster is quite clumpy, and thus a centroid cannot be uniquely defined. That said, at a distance of 
100" from the brightest cluster member, we find a tangential shear of about 0.2, in excellent agreement with the weak 
lensing study of L06. 

We have produced a simple shear field as described in the previous section and used this field to reconstruct the density 
field via the Seitz & Schneider (1995) smoothed finite-inversion technique. Our shear field and reconstructed k field may 
be seen in Figure [3] 

With so much attention given to the ability of strong lensing to pick out substructure in clusters, it is gratifying to 
note that weak lensing alone can identify large concentrations of galaxies. A contour plot of the projected surface mass 
density, k, is shown in Figure [4] overlaid on an image of the cluster itself. 
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Fig. 3. — The weak-lensing shear field along with a non-parametric reconstrucion for ACS images of A1689. Note that 
the overall normalization of the reconstruction is subject to uncertainty due to the mass-sheet degeneracy in a finite field. 
North points to the lower right of the image, and the width of each bin corresponds to approximately 45/i _1 kpc at the 
distance of the cluster. 



For this cluster, we find X cr it ~ 5.7 x 10 15 /iM Q Mpc -2 (taking z = 0.18), which allows us to estimate the true surface 
mass density up to the mass-sheet degeneracy. However, since we have not independently measured the shear field far 
from the center of the cluster, we rely on the estimates of others. L06, for example, found a mean k of 0.48 within 1' of 
the center of the cluster. We have thus set our mass sheet at a slightly lower threshold, so that the mean over the entire 
chip is 0.4. 

We note that the density profile is quite shallow, however. We find that a power law profile is well fit if a = —0.19. 
Over a similar range, L06 fit to a profile with a ~ —0.6. L06, however, base their fit on extrapolation from the slope of 
the shear outside the ACS image. Inspection of Fig. 11 in their work suggests that the circularly averaged profile within 
the central arcminute may be much shallower. 

4.2. Strong+Weak Lensing Mass Reconstruction (Non-Parametric) 

We have combined strong lensing information from L06 with our shear measurements to generate a combined convergence 
map, following the method of Bradac et al. (2005a). A convergence map and a contour plot of the convergence found by 
this method can be seen in Figures [5] and [6] 

Figure [6] clearly shows an elongation in the direction of the secondary dark matter clump described in L06, which is 
not seen in Figure g) Figure |7] compares the circularly averaged convergence as a function of distance from the center of 
the cluster for the weak and strong+weak mass reconstructions. The strong+weak profile has a much steeper slope than 
the weak profile. Indeed, we find a — —0.29 for the profile shown in Figure [7] 

This difference in slope could result from the fact that our masking scheme results in an under density of background 
sources in the central region of the cluster. This shortage of data points results in a lower shear signal in the central region, 
thus lowering the computed value of k. The under density of sources in the central region is found to be an important 



factor in the errors associated with our parametric flexion reconstruction, and is discussed in more detail in § 4.3 
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4.3. Flexion Mass Reconstruction (Parametric) 

A parametric mass reconstruction was generated using (first-) flexion data alone. In dense systems like clusters, we 
have found that the HOLICs technique is less susceptible than shapelets to contamination by light from the extended 
wings of lens galaxies and other neighboring sources. Flexion probes the higher order shape moments, so even a small 
contamination near the edges of postage stamps can cause a signicant spurious flexion signal. We thus used our HOLICs 
measurements exclusively to estimate the flexion signal. 

The reconstructed convergence is plotted in Figure [8 and a contour plot is shown in Figure [9] The reconstruction 
shows significant substructure, which appears to be well correlated with small clumps of galaxies outside the center of 
the cluster. However, there is a rather worrying under density seen in the center of the cluster. In order to asses the 
significance of this under density, it is necessary to quantify the errors associated with the convergence map. 

Figure [10] shows the approximate errors in the flexion reconstruction. These errors were computed as follows: each 
flexion data point was rotated by an angle drawn from a uniform random distribution in the range [0, 2tt]. Parametric 
reconstructions were generated using this randomized data, and this procedure was carried out for 1000 randomizations. 
The errors presented in the figure are the rms values found for each bin. 

Clearly, the under density seen in the center of the convergence map in Figure [8] should be considered as resulting from 
noise, rather than a real feature of the cluster, as the central region appears to be entirely noise-dominated. This is most 
likely due to the fact that there are fewer background sources found in this region. Figure [TT] shows the average number 
density of background sources plotted as a function of the radius over which this number density is averaged. 

The shortage of sources in the central region results from our masking scheme, and affects our flexion reconstruction in 
two ways. Firstly, it means that a foreground galaxy in the center of the cluster is less likely to have the required 3 nearby 
background sources, and thus fewer of the central cluster members will be included in the analysis. Secondly, those that 
are included will generally be fit using fewer data points than the outlying cluster members, and thus these fits will have 
larger associated errors. 

Thus, the under density seen in the center of the flexion reconstruction should not be believed. However, in the outlying 
regions of the image, the noise is seen to drop significantly, and the substructures seen in these regions appear to be real 
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features. 

4.4. Galaxy- Galaxy Flexion Signal 
In addition to a large-scale map of the clus ter, we have also generated a composite circular profile of the cluster member 



galaxies via flexion measurements. Figure 12 shows the galaxy-galaxy (first-) flexion averaged over all background- 
foreground pairs meeting our selection criteria. 

We find that the mean cluster member can be fit well by an isothermal sphere with o~ v = 295 ± 40 km/ s. Given 
the relatively large uncertainty in velocity, and the relatively narrow scatter in the magnitude of the cluster members 
(<r r ~ 1), we are unable to effectively split the cluster members into subgroups. However, we can compare this result to 
that expected from the Faber-Jackon relation. 

The mean absolute magnitude of the sample is approximately -21.5 in the R-band. Taking the canonical Faber-Jackson 
relation: 

L ^ 25 



a v = 220 km/s 1-1 , (28) 

we find an expected velocity dispersion of approximately 390 km/s, in accord with our measurements. 

L06 (Table 3) compute a best estimate RMS for several member galaxies (~ 200km /s) and the Brightest Cluster Galaxy 
(~ 500fcm/s). Our mean flexion estimate falls squarely in the middle of this distribution. 

We also note that the mean velocity dispersion found when computing the parametric flexion reconstruction (321fcm/s) 
is consistent with that derived from the galaxy-galaxy lensing study, within the error bars of the latter. 

5. SUMMARY 

We have used a shapelet-based shear measurement technique to create a non-parametric mass reconstruction of the 
galaxy cluster Abell 1689. Using only weak lensing data, we found significant ellipticity and substructure in the cluster. 
Combining the weak lensing data with strong lensing data from L06 improved the resolution of the mass map, and 
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Fig. 6. — A contour plot of the convergence, k, determined using the shapelet measurements of the ellipticity of background 
sources in the field combined with strong lcnsing data from L06. 
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Fig. 7. — The circularly averaged convergence as a function of distance from the center of the cluster from weak 
measurements alone (dashed) and from the combined weak+strong calculation (solid line). Clearly, the weak+strong 
profile shows a much steeper slope than the weak profile. However, there is very good agreement between the two profiles 
at large radii. The error bars shown in these plots are taken from y/N statistics, thus are only an approximation. 

increased the slope of the cuspy central density profile. The combined analysis also identified a secondary dark matter 
clump found by L06. Using an entirely new and independent flexion analysis, we were able to verify the position of 
this clump, and other substructure, via a parametric reconstruction of the cluster mass, modeling each of the cluster 
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Fig. 8. — The binned convergence found using a parametric model for the cluster galaxies from flexion data. The width of 
each bin is approximately 22.5/i _1 kpc in the cluster plane. Note that a flexion-only reconstruction is sensitive to localized 
substructure but not smooth gradients in the mass distribution. 

members as a singular isothermal sphere. The substructure observed in the flexion reconstruction is well-correlated with 
the locations of groups of cluster galaxies. 

We have also used flexion data to probe the halos of individual cluster galaxies. Using a similar parametric recon- 
struction, we measure a highly significant (~ 13c) galaxy-galaxy flexion signal. In agreement with previous, non-flexion 
measurements, we find a mean velocity dispersion of 295 ± 40 km/s for the cluster galaxies. 

The authors would like to thank David Bacon, Sanghamitra Deb, John Parejko, Lindsay King, and Marceau Limousin 
for useful discussions. This work was supported by NASA ATP NNG05GF61G and HST Archival grant 10658.01-A. AL 
is supported by a BP/PPARC Dorothy Hodgkin Postgraduate Award. 

APPENDIX 

INDUCED FLEXION DUE TO THE POINT SPREAD FUNCTION 
From equations 43-45 in GL07, we have the relations: 

Qijk = Qijk + Pijkt 

Qijki=Qi%i + dP ijk i, (Al) 

where 

dPijki oc P ijkl cx ap SF (A2) 
anC ^ (o) 

Qijkl ^ a source- (A3) 
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Fig. 9. — A contour plot of the convergence calculated using a parametric model for the cluster galaxies from flexion 
measurements. This reconstruction shows significant substructure corresponding to the locations of small clumps of cluster 
galaxies. 



Now 

T 



Qijk 
Qijkl 

Qijkl + dP ijkl 

n (0) n (0) u av 

^ijk ^ijkl *ijk U-Lijkl 



~ ^ 4 <OU : C \ +?PSF 4 4 ■ (A4) 

a source ' a PSF a source ' a PSF 

The former term in the above expression will dominate for a source > apsF, and the coefficient of the J 7 ^ term will 
approach unity as a source becomes large. In the case of small a source , the term in Tpsf will become important, provided 
that FpsF itself is significant. 
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Fig. 10. — The figure shows a contour plot of the approximate errors in the parametric flexion mass reconstruction. It is 
evident from this plot that the reconstruction is entirely dominated by noise in the center of the image. 
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Fig. 11. — The figure shows the number density of background sources used in our flexion study (averaged over a circular 
region) plotted as a function of distance from the center of the image. The dotted line shows the number average number 
density over the entire image. There is an apparent downturn in the number density at large r. This downturn occurs 
simply because we reach the edge of the field. 
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